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Abstract 


A general-purpose digital computer program (named ELAS) for the in-core 
solution of linear equilibrium problems of structural mechanics is described for 
potential and actual users in Volume I of this report, and is documented in 
Volume II. The program requires minimum input for the description of the 
problem. The solution is obtained by means of the displacement method and 
the finite element technique. Almost any geometry and structure may be handled 
because of the availability of lineal, triangular, quadrilateral, tetrahedral, hexahe- 
dral, conical, triangular torus, and quadrilateral torus elements. The assumption 
of piecewise linear deflection distribution insures monotonic convergence of the 
deflections from the stiffer side with decreasing mesh size. The stresses are pro- 
vided by the best-fit strain tensors in the least-squares sense at the mesh points 
where the deflections are given. The selection of local coordinate systems when- 
ever necessary is automatic. The core memory is efficiently used by means of 
dynamic memory allocation, an optional mesh-point relabelling scheme, imposi- 
tion of the boundary conditions during the assembly time, and the straight-line 
storage of the rows of the stiffness matrix within variable bandwidth and the 
main diagonal. The number of unsuppressed degrees of freedom that can be 
handled in a given problem is 500 to 600 for a typical structure, but might far 
exceed these average values for special types of problems; the execution time 
of such problems is about four minutes in 32K IBM 7094 Model I machines. The 
program is written in FORTRAN II language. The source deck consists of about 
8000 cards and the object deck contains about 1400 binary cards. The physical 
program (standard ELAS) is available from COSMIC, the agency for the dis- 
tribution of NASA computer programs. 


vi 
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I. Introduction 


Volume I, Users Manual , of this report gives a general 
description of ELAS,* a general-purpose digital computer 
program for the in-core solution of linear equilibrium 
problems of structural mechanics, and contains the infor- 
mation necessary for input preparation, arrangement of 
the physical program, and interpretation of output and 
error messages. 

Volume II, Documentation of the Program , is published 
in two parts: the present volume— the basic Volume II— 
which gives the theoretical background of the program 
and contains tables and figures describing the COMMON 
variables, their meanings, and their arrangement in 
COMMON; and an Addendum to Volume II, which 


*First two syllables of the word Elasticity. 


contains program descriptions, flowcharts, and source 
program listings for all program elements of ELAS/Level 
3. (The original version of the ELAS program made 
available from COSMIC** in April 1968 is designated 
ELAS/Level 0. Subsequent program corrections made in 
January 1969, March 1969, and May 1969 updated the 
program to ELAS/Level 1, ELAS/Level 2, and ELAS/ 
Level 3, respectively.) 

In addition to the list of references cited in the text, 
a list of documented works related with the development 
of the ELAS program is given in the bibliography. A 
corrigenda for Volume I is given in the Appendix. 

**Computer Software Management and Information Center, Com- 
puter Center, University of Georgia, Athens, Georgia, 30601, 
telephone 404-542-3265. 
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II. Theoretical Background 


This section summarizes the mathematical formulation, 
the numerical method of solution, and the design features 
of the program. 

A. Mathematical Formulation 

Let V denote the material volume of the structure 
within the closed boundary S. Let x a , a = 1, 2,3, denote 
a fixed right-handed Cartesian coordinate system. The 
Greek subscripts always refer to these axes; therefore, a a p 
is the stress tensor described in such a coordinate system. 
Let u a denote the displacement vector; p a the body force; 
m the unit mass; double dots above, the second time de- 
rivative; and comma in the subscript the partial differen- 
tiation with respect to the space variable represented or 
implied by the subscript following the comma. If it is 
assumed that repeated subscripts imply summation over 
the range, the equilibrium of any particle within V may 
be expressed as 

o- 0 a >/3 + Pa = ntiia (I) 

In the equilibrium problems, the loading is such that 

Ua - 0 (2) 


Therefore, substitution of u a from Eq. (2) into Eq. (I) 
yields 

+ pa = 0 in V (3) 

Let S' denote the portion of S where the tractions are 
prescribed. The equilibrium condition on S' is 

(Tap Up + Pa = 0 (4) 

where n$ is the unit normal vector and p a is the prescribed 
traction. The stress-strain relationship of the material is 

DygyS (^yg € yg) (5) 

where e^ d is the prescribed strain tensor, € y6 is the strain 
tensor, and D a pys is the material matrix, which is positive- 
definite and symmetrical, so that 

Dap yd — Dy$ap — DpayS ~ Dapdy 

The strain displacement relationships are 

€«3 = | + ( 6 ) 
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Finally, the displacement boundary conditions may be 
stated as 

«* = «£ onS" (7) 

where u° a denotes the prescribed displacements on S". 
It should be noted that 

S' + S" - S (8) 

In an equilibrium problem, usually V, S', S", p a , u° a , p a> 
Dapys, np, and e° 6 are given and w«, € a p, and <r a p are 
requested. 

Equations (3) through (8) constitute the differential 
equation formulation of the equilibrium problem in three- 
dimensional continuum. A finite difference solution based 
on this formulation may be set up as follows. A regular 
mesh is placed in V such that S' and S" are determined 
by the mesh points. If S is not defined by coordinate 
surfaces, such representation of S' and S" is only approx- 
imate. The displacements u a at the mesh points in V and 
on S are taken as the primary unknowns. The prescribed 
u° a displacements are assigned to the mesh points of S". 
With the use of Eq. (6), e a p at the mesh points in V and 
on S are approximated by the first differences of u a and 
u°. Then, by the use of Eq. (5), the values of a a p are 
expressed at the mesh points. Finally, depending upon the 
mesh point in V or on S', Eq. (3) or Eq. (4), respectively, 
is used to write the difference equations for the unknown 
displacements. After the unknown displacements from 
these equations have been computed, the strains and the 
stresses may be computed from the finite difference ap- 
proximations of Eq. (6) and Eq. (5). Such a solution 
method has the following drawbacks: 

(1) To minimize the truncation errors, a regular mesh 
in V is used; however, this causes approximate 
representation of boundary S and, therefore, in- 
creases the truncation errors in the finite difference 
approximations of Eqs. (4) and (7). Since the errors 
in the finite difference approximations of Eqs. (4) 
and (7) dominate in the solution more (Ref. 1) than 
the errors in the finite difference approximation of 
Eq. (3), either an irregular mesh in V may be 
considered to represent S more accurately, or 
higher-order formulas for the boundary conditions 
are used, although neither scheme is desirable in 
a general-purpose digital computer program. 

(2) Because of the symmetry and the positive-definite- 
ness of D a py8, the formulation given by Eqs. (3) 
through (8) is self-adjoint and positive-definite. 
However, the coefficient matrix of the unknown 


displacements in the finite difference equilibrium 
equations is, in general, neither symmetric nor 
positive-definite. The loss of the two desirable 
qualities of the problem in the numerical formula- 
tion increases the storage requirements and solution 
time. Because of these setbacks, the mathematical 
formulation given by Eqs. (3) through (8) is modi- 
fied slightly as explained in the following paragraph. 

Consider the quantity ir defined as 
7 r = 2" &ctp dV ^ U® Pa (TV J’ tla Pa dS 

(9) 

where dV is the volume element, dS the area element, 
and the other symbols are as previously defined. Consider 
the displacement fields satisfying Eq. (7). For each such 
displacement field, by means of Eqs. (9), (6), and (5), a 
scalar tt may be computed. It can be shown that, for 
sufficiently smooth displacement fields satisfying Eq. (7), 
the stationary point of tt, i.e., the point for which 

S 7T — — 0 (10) 

also satisfies Eqs. (3) and (4). In fact, by the methods of 
calculus of variations, Eq. (10) yields Eq. (3) as the Euler 
differential equation, and Eq. (4) as the additional bound- 
ary condition. Therefore, Eq. (10) is an equivalent state- 
ment of Eqs. (3) and (4). The quantity tt is known as the 
“total potential energy” of the system. Thus, the formula- 
tion given by Eq. (10) reduces the problem to that of 
locating the stationary point of the total potential energy 
functional. How the numerical solution is set up from 
this formulation (which is sometimes referred to as the 
extremum formulation of the problem), and its advantages, 
are discussed in the next subsection. 

B. Numerical Solution Method Based on the 

Extremum Formulation 

A random mesh is placed in V such that the mesh 
elements are line segments, triangles, quadrilaterals, tetra- 
hedrons, hexahedrons, conical segments, or triangular or 
quadrilateral tori. Some of the mesh elements are shown 
in Fig. III-l (Vol. I). The types of mesh elements that may 
be used in different structures are given in Table III-2 
(Vol. I). The randomness of the mesh enables the selection 
of mesh points that are exactly on the boundary S. For 
clarity, the mesh points are labelled sequentially, with 
integer numbers starting from 1. If there are £ number of 
mesh points, there are s! different types of possible label- 
ling. In the discussion that follows it will become obvious 
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that some of these labelling systems are more desirable 
than others. 

It is assumed that one of the possible s! systems is 
selected. Next, the mesh elements are labelled sequen- 
tially with integers. If there are m number of mesh ele- 
ments, there are ml number of different labelling systems. 
It is supposed that one of the possible m! systems is 
selected. In what follows, superscript m indicates the 
element label, and subscripts t or s indicate the mesh 
point label. To solve the equilibrium problem formulated 
in Section II-A numerically, instead of computing u a at 
every point of V and S, an attempt is made to find, at the 
mesh points, certain related quantities that define the 
distorted configuration of the structure in the same way 
as u a . These quantities are called deflections, which are 
displacements /rotations at the mesh points. Given a mesh 
point, the total number of independent deflection com- 
ponents is the number of degrees of freedom of that mesh 
point. In Table III-l (Vol. I), the deflection components 
at a mesh point of different structures are shown as re- 
ferred to an overall coordinate system. Let k denote the 
number of degrees of freedom at a mesh point of a struc- 
ture. The value of k for different structures is given in 
the last column of Table III-l (Vol. I). It will be assumed 
that the deflection components at a mesh point are ordered 
as shown in the table. The subscripts k and l will be used 
to indicate the sequence number implied by this ordering. 
If a prime appears on k or Z, this implies that a local co- 
ordinate system is used in defining the degree-of -freedom 
directions. Let q k S denote the kth deflection component 
at mesh point s. A mesh element may be defined by the 
mesh points that are coincident with its vertices. For 
clarity in referencing, the convention of Table III-5 
(Vol. I) is adopted in ordering the vertices of mesh ele- 
ments. The type numbers shown in this table refer to the 
numbers shown on the shaded squares of Table III-2 
(Vol. I). Subscripts g and h will be exclusively used to 
denote the sequence number of a vertex in the g number 
of vertices of an element. 

With the preceding definitions, the method used to 
obtain the stationary point of the total potential energy 
functional may now be explained. This is the classical 
Ritz procedure (Ref. 2), where the undetermined param- 
eters of the problem are the unknown components of q^s 
deflections. Equation (9) is first written as 


where the repeated index implies the summation over the 
range. Next, an attempt is made to select a family of 
displacement fields that are sufficiently smooth, but 
otherwise arbitrary, ignoring for the time being the 
essential boundary conditions of Eq. (7). A piecewise 
linear displacement field is acceptable in this sense (Ref. 3). 
Of course there are other piecewise continuous fields that 
are acceptably smooth. However, to simplify the under- 
standing of the procedure, it is assumed that the displace- 
ment fields are piecewise linear. Such a field may be 
described mathematically for the mth element in terms 
of the deflections of its vertices as 

'V' 

~ B a- fi- I'hQi'i Put <ht V + rigid body movement 

( 12 ) 

where the primes indicate the local coordinate system 
of the element. The coefficients fj^ t constitute a binary 
array such that, for a given m and h , it is zero throughout 
the range of £, but 1 at the value of t corresponding to 
the hth vertex of the mth element. In fact, q u is the 
list of deflection components pertaining to the vertices 
of the mth element. The matrix Qu i in Eq. (12) is the 
coordinate transformation matrix, where, for fixed l\ it 
represents the direction cosines of the local axis related 
with V degree-of-freedom direction in the coordinate sys- 
tem associated with l. The space variable is the dis- 
tance measured along the /?' th local axis. The coefficients 
Bct'p'i'h may be computed from the local coordinates of 
the vertices of the mth element. In Table III-3 (Vol. I), 
for different mesh elements, the orientation of the local 
coordinate system relative to the overall coordinate sys- 
tem is given. It should be noted that Eq. (12) is an 
approximation of the true displacements in the mth ele- 
ment, even if the exact values of deflection components 
qu are known. However, it may be shown that the error 
decreases with decreasing mesh size. With the use of u a > 
from Eq. (12) in Eq. (6), the strains in the mth element, 
as referred to the local coordinate system of the mth 
element, may be expressed as 

<s/ 

= Boi'p'l'h Ql l ! x *httylt (13) 

Let Da'P'y 8' denote the material constants of the mth 
element as referred to the local coordinate system of the 
mth element. If 


" f j Ua Pa dV m f Ha pa dS m 

" J V™ jy« JS'n 

( 11 ) 


K| 


kglh 




B». 


r'fc ' 


,D 


g'r'a'P 


' Ba’ 


P’l' 


hdV ) Qi • i 

(14) 
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and 


Pftg — Qk'k^J &a' fi'.k' ffXp' Pa' 


are defined, it of Eq, (11) may be expressed as 

77 = *2 9*3 V'Ts^Tgih Put 9 It ~~ Qks (p^s^fg Qks) ( 16 ) 

where Q ks denotes the prescribed concentrated loads at 
the mesh points. The deflection components q fcs (or qu) 
should be such that, on S', they satisfy the essential bound- 
ary conditions of Eq. (7). Let dj denote the portion of qu, 
which is unknown. The essential boundary conditions may 
be expressed as 

qit — eitjdj 4* e° lt (17) 

where the coefficients eu$ and e° lt are quantities that may 
easily be determined from Eq. (7). For example, if there 
are no prescribed deflections in the problem, ei t j is a 
binary array containing only one 1 in the whole range of 
/ for a given It, and e° lt is zero throughout. Actually, the 
dj are the true undetermined parameters of the problem. 
If qu is substituted from Eq. (17) into Eq. (16), the values 
of dj corresponding to the stationary point of it may be 
obtained from the set of linear equations 

•= 0 (18) 

since 

§7T = 7 T } dj $dj (19) 

The equations given by Eq. (18) may be rewritten as 

&Cn dj = SPi (20) 

where 

9iij = &ksi (Ajs^kglh ^ht e itj ( 21 ) 

and 

“ &ksi {^gs^lzg “b Qks) ~~ &ksi ^gs^kglh P Tt e °lt (22) 

The coefficients Kf glh constitute the element stiffness 
matrix of the rath element and P™ g is called the rath ele- 
ment load vector. In Eq. (20), the coefficient matrix QCa 
is the stiffness matrix associated with the directions of dj 
deflections, and the right-hand-side vector 0i lists the 
loads in these directions. Equations (21) and (22) indicate 
how the coefficient matrix and the right-hand-side vector 
of the governing equations can be systematically gener- 
ated from the element stiffness matrices and load vectors. 


dV Ba' ft'k' gXp' p U ' dS J (15) 


The operation implied by Eqs. (21) and (22) is referred 
to as the assembling of the elemental matrices. Because 
of the positive-definiteness and the symmetry of D a p y d, 
Eq. (21) shows that 

Qlij — Qiji (23) 

and QUj is also positive-definite. Once the unknown de- 
flections dj are solved from Eq. (20), the complete deflec- 
tions qu are obtained by substituting dj into Eq. (17). 
After deflections q ks have been computed, the strains and 
the stresses at the mesh points may be computed as 
described in Ref. 4. 

The method of solution described in the preceding has 
the following advantages: 

(1) Since the mesh is random, the boundaries S' and S" 
may be closely approximated, and thus minimize 
truncation errors. 

(2) Any a priori knowledge about the variation of u a in 
V or on S may be used advantageously by varying 
the mesh size accordingly to minimize the trunca- 
tion errors. 

(3) The self-adjoint character, as well as the positive- 
definiteness of the problem, is preserved, since 9ijj 
is always symmetric and positive definite. 

C. The Program 

I. Criteria for Storage Allocation . AH the input data 
previously mentioned (see Fig. IV-1, Vol. I) are stored 
permanently in COMMON after their validity is checked. 
No fixed-length block is assigned to these diversified data. 
The data are compactly stored as a string of variable 
length. This enables the program to compete in obtaining 
priorities efficiently in a multiprogramming environment. 
After reading the control card, the program determines 
the pointers of each of the data blocks relative to the 
beginning of COMMON. When the other input data 
become available, they are placed in the proper place in 
COMMON by the pointers. Although the locations of the 
data blocks vary from one job to another, the locations 
of the pointers and the control information provided by 
the control card have fixed locations at the beginning of 
COMMON. The remainder of the core is assigned for 


6 


JPL TECHNICAL REPORT 32-1240 



the program instructions and temporary storage for the 
coefficient matrix and the right-hand-side vector of the 
governing equations. The program consists of four links, 
and since all the program instructions are not required 
simultaneously, only the instructions for each link in turn 
as needed are retained in the core. 

A sketch of the governing equations in Eq. (20) is given 
in Fig. II-l (Vol. I). Since the coefficient matrix is sym- 
metric, the program allows storage only for the shaded 
area shown In the figure. From Eq. (20) it may be ob- 
served that, for a fixed /, 9C%$ represents a vector listing 
the forces that may develop when a unit deflection is 
applied in the fixed / direction, while keeping all the other 
degree-of-freedom directions with zero deflections. The 
nonzero entries of this vector coincide with the deflec- 
tion directions of only those mesh points that are con- 
nected with the disturbed mesh point by means of mesh 
elements and deflection boundary conditions. This shows 
that the Qii) matrix is sparse and usually has a large zero 
area in the upper right-hand comer. 


Before generating the coefficient matrix and the right- 
hand side of Eq. (20), the program computes a pointer 
for each of the rows of the coefficient matrix, and a pointer 
for the right-hand-side vector so that the coefficients 
shown 'in the shaded area in Fig. II-l (Vol. I) can be 
stored compactly in COMMON as a string. Actually, the 
pointers of the rows are the addresses of the words imme- 
diately preceding the diagonal elements. As discussed in 
Ref. 5, by the proper ordering of dj unknowns in Eq. (20), 
the zero area may be increased in the upper right-hand 
comer of Qii If the user chooses to assign zero into the 
ISHUF field of the control card, the unknowns are 
ordered as implied by the mesh-point labels; e.g., the 
unknown deflection components of the first mesh point 
are placed first, those of the second mesh point are placed 
second, etc. If the user assigns ISHUF = 1 or 2, the pro- 
gram first tries to find a better labelling system with the 
method given in Ref. 5, and uses these new mesh-point 
labels in ordering the unknowns dj. For example, if mesh 
point with label 25 is the first mesh point in the new 
labelling system, the unknowns of this mesh point are 
listed first in dj. If the user assigns ISHUF = 3, the better 
labelling system is required by the program from input 
data cards (number 17 in Fig. IV-1, Vol. I). The method 
for relabelling described in Ref. 5 requires the genera- 
tion of the mesh-point connectivity matrix N s t, which is 
a binary matrix. If mesh point s is connected to mesh 
point t by a mesh element or by a deflection boundary 
condition, N 8t = N ts — 1; otherwise, N st = N ts = 0. It is 
always assumed that a mesh point is connected to itself. 



Fig. 11-1. Definition of Eiy augmented matrix 


If point s is completely constrained by the deflection 
boundary conditions, N st = N ts = 0 for all t, except t = s. 
The program generates N st from the information pro- 
vided by the element description data and deflection 
boundary conditions. The connectivity matrix N st is 
always generated, since it is also used in determining 
the pointers of the rows of Qiiy. 

2 . Method of Assembly . To obtain the coefficient matrix 
and the right-hand-side vector of the governing equa- 
tions, the mesh elements are processed, one at a time, 
first to obtain the element stiffness matrix and the ele- 
ment load vector for each, and then to assemble these 
according to Eqs. (21) and (22) and the allocated storage. 
Let A r denote the vector in COMMON and R* denote 
the pointer of the tth equation in Eq. (20). Let us assume 
that the right-hand-side vector is stored after the coeffi- 
cient matrix. Let E ksi > (or Eity) denote the augmented 
matrix composed of e kS i (or e U y) and ef iS (ore° lt ), as shown 
in Fig. II-l. Let a prime on the index indicate that the 
range of unprimed index is increased by 1, and let an 
underlined index indicate the largest value within the 
range. With this notation, the assembly procedure may be 
summarized as 


A r — y i'Eksi'Qks + 

8i’j‘E ksi , p,™K™ sll h (25) 

where 

8 l.j, = 1 if i' s=/', f and r — Rj< + f — i'\ 

8l,y = — 1 if f >i, i, andr = Rj -f *' > (26) 

8? (/ , = 0 for all other possibilities / 

and 

y\, = 1 if i' and r = Ri + i' f 

> (27) 

y\, = 0 for all other possibilities ) 
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In the program, only the nonzero E k8i * constants are 
computed by the deflection boundary condition input 
units and the connectivity matrix. For each ks the non- 
zero entries of E ksi r are stored with their values and V 
indices. The values and the indices of nonzero Q ks entries 
are directly provided by the concentrated load input 
cards. The binary coefficients 8 r irjf and y\, are not stored, 
but determined from Eqs. (26) and (27). If m and g are 
given, the s value of the nonzero entry of /a™ is obtained 
from the element description data. Let el and i' a denote 
the nonzero values and corresponding indices in E k8i * 
for a given ks, and let a denote the maximum value of a, 
so that l^a^a ( b and b are alternate symbols). Let d 
denote the number of concentrated load input units. This 
notation is used in the flow diagram corresponding to 
Eq, (25) given in Fig. II-2. In the ELAS program, the 
summations implied by the first term in the right-hand 
side of Eq. (25) are implemented in Link 1 and the re- 
mainder in Link 2. 

3. Method of Solution of the Governing Equations. 

Since QUj is a symmetric and positive-definite and band- 
limited matrix for the solution of Eq. (20), the Cholesky 
algorithm may be applied. In this method, the decom- 
posed stiffness matrix Bij from 9i%j is first computed as 

Bij * Bj'j ~ 9(ij (28) 

where the range of f equals that of i. Then, from 

Bijd'j — (29) 

with a forward sweep, the auxiliary unknowns d' can be 
solved. Finally, 

Bij di = d'j (30) 

yields the unknowns with a backward sweep. In Fig. II-l 
(Vol. I) the border of the zero area in the upper right- 
hand corner is not always defined by the last nonzero 
coefficient in each equation. This is because the shaded 
areas of Qlij and Bij are identical only when the border is 
selected as shown in Fig. II-l (Vol. I). In the ELAS pro- 
gram, 0Uj coefficients in the shaded area of the figure are 
first modified to those of the coefficients of B i; *, then ffi 
constants are changed to d\ , and finally, di values are con- 
verted to the numerical values of the dj unknowns. Then, 
from Eq. (17), q u is computed on the same area as dj. 
These operations are carried out in Link 3 of ELAS. 

4. Computation of Stresses. The computation of stresses 
in displacement methods poses a harder problem in struc- 
tures of two- or three-dimensional continuum than that in 
truss and frame structures, which truly have a finite num- 



Fig. 11-2. Flow diagram corresponding to the summations 
implied by Eq. (25) 
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ber of deflection components for the determination of 
their distorted configuration. The difficulty arises from 
the fact that the structures of two- or three-dimensional 
continuum actually have infinitely many deflection com- 
ponents, and the relations of the type of Eq. (12) are only 
approximate. 

After computing the deflections as mesh functions, the 
problem of stress computation with acceptable accuracy 
in reasonable machine times still remains. Experience has 
shown that the use of Eq, (13) and then Eq. (6) presents 
the following drawbacks: (1) the exact location of the point 
for which the stresses are computed is not known, and 
(2) the computed stresses may be largely different from 
the actual stresses. Despite these setbacks, stress compu- 
tation of this type is being widely used because it is modu- 
lar in elements, just as is the generation of the governing 
equations in Eq. (20), a feature that facilitates automa- 
tion. In the ELAS program, the best-fit stress computation 
method of Ref. 4 is used for structures of continuum. This 
method is just as easy to automate and has the following 
advantages: (1) stresses are computed at the points where 
the deflections are obtained, (2) the accuracy in stresses 
is comparable with that of deflections, and (3) the stress 
boundary conditions of Eq. (4) may be satisfied during 
the computation of the stresses of boundary points. This 
scheme was initially devised for triangular finite elements 
(Ref. 6). 

In the following paragraphs the stress computation 
in structures of two- and three-dimensional continua is 
explained. The computation of stresses in structures com- 
posed of elements of one-dimensional continuum is per- 
formed by multiplying element stiffness matrices with 
computed deflections. 

Mesh Line Set. Suppose that the deflections at the mesh 
points of a structure of three-dimensional continuum are 
available and that the stresses at mesh point s are re- 
quested. The question of how much deflection data should 
be included in the computation is of practical importance 
because the computation time rapidly increases with this 
quantity. Experience with the method of computing 
stresses in the element indicates that deflections of the 
set of elements meeting at mesh point s are sufficient for 
the computation of its stresses with acceptable accuracy. 
The mesh points of the element set are called “mesh-point 
set” and the mesh lines meeting at the common mesh 
point s are called “mesh-line set.” The scheme adopted in 
ELAS is modular in the mesh-line set— the next-best unit 
after elements. The stress computation at a mesh point 
starts with the determination of the element set, and 


consequently, the mesh-line set associated with this 
mesh point. Then, if this mesh point is on the bound- 
ary, the average boundary surface area associated with 
this node and the direction cosines of the outer normal 
are computed. 

Selection of Local Coordinate Systems at the Mesh 
Points. In a given problem, it is desirable to have one 
fixed, right-handed coordinate system to express the 
stresses. However, this is not practical for structures com- 
posed of anisotropic material, at the boundary points 
where the outer normal is not coincident with the coordi- 
nate lines, and for shell structures. The following method 
is adopted in the ELAS program for the selection of local 
coordinate systems at the mesh points. 

At an internal node, the local axes may be taken as the 
material axes unless the material is isotropic, in which 
ease they should be taken (1) parallel to the overall 
coordinate system in plates and three-dimensional solids, 
and (2) as the principal curvature directions and the nor- 
mal of the middle surface, or any other suitable system 
that the user inputs in shells. At a boundary node, the 
first local axis may be coincident with the outer normal, 
and the directions of the remaining local axes may be 
determined so that (1) the local third axis becomes the 
middle surface normal in plates and shells, and (2) the 
direction defined by the cross-product of outer normal 
with the overall axis, which makes the largest angle with 
the outer normal, then becomes the second local axis in 
three-dimensional solid structures. 

Stress Computation at an Internal Mesh Point . Let / be 
the label of a mesh line in the mesh-line set, with A py 
the position vector and A the displacement vector of 
the far end of the mesh line relative to the mesh point 
where the local coordinates are defined. As the first 
approximation of the strain along the ;th mesh line, the 
following may be written: 

€ = (31) 

a Pfi' Apr 

The same strain may be obtained from the strain tensor 
of mesh point s as 

e= Apg.Ap^lg: (32) 

Apr A 

Equating Eq, (31) to Eq. (32) and cancelling the denomi- 
nators results in the following expression: 

(Ap<*' A p&')j €at'fS' = (Apy' A Uy') j (33) 
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The number of equations in Eq. (33) is equal to the 
range /. Usually the range of / is greater than the num- 
ber of independent components of the strain tensor. (If 
not, the mesh may be readjusted by repeating the deflec- 
tion computation.) Therefore, in Eq. (33), there are more 
equations than the unknown strain components. Such a 
set may be solved by least squares, first by multiplying 
both sides with the transpose of the coefficient matrix, 
then by inverting the new coefficient matrix. This leads to 

€cc'P' = [(A p & ' Ap$')i(Ap5' Apv')i] -1 (Apfi/ Apv'); (Apy' A Uy')j 

(34) 

where the range of i is equal to that of f. The stresses at 
the mesh point may be obtained by substituting from 
Eq. (34) into Eq. (5). If the problem is a plane-strain 
problem, one should first impose 

£3 '3' ~ € 1'3 # — € 3 / ; l ' = 62' 3' ^ €3*2' ~ 0 (35) 

on Eq. (33). For plane-stress problems, D a ^ yd' in Eq. (5) 
should be modified to guarantee 

Ci'3' = G&'I’ “ € 2 '$' = € 3 * 2 ' — <JV3' == 0 (36) 

For the bending of plates and shells, e ar ^ should be inter- 
preted as curvature changes and, in Eq. (33), (ApyvA%^) ; - 
should be taken as the projection of the rotations vector 
of the far-end mesh point relative to the current mesh 
point in the /th mesh line on X A p direction, where % is 
the unit vector of the third local axis. Also, the conditions 
in Eq. (36) should be imposed on and Eq. (5) 

should be replaced by 

t 3 

My'd' 22 Dy'5'a'fi' €«'/}' (37) 


where M 7 ' 5 ' denotes the bending moments and t is the 
thickness. The membrane case of shells is identical with 
the plane-stress case, provided that Eq. (5) is replaced 
with 

Ny'd' ~ t jDy'g' a'fi' (38) 

where Ny'S' denotes the membrane forces. 

Stress Computation at a Boundary Mesh Point. The 
procedure for stress computation at a boundary mesh 
point is basically the same as the computation at an 
internal mesh point. Here, the stress boundary conditions, 
expressed in terms of the strains, are included in Eq. (33) 
before the application of the least-squares scheme for 
their solution. The stress boundary conditions may be 
written as 

j ^ Pa' /'xn 

n Da'l 3 ' 7 ' 6 ' € 7 ' 6 ' — r> (3 y ) 

where represents the prescribed boundary stresses. 
If at the boundary mesh point the deflections, in place 
of the stresses, are prescribed, R a ' reaction forces may 
be found from the equilibrium equations of the boundary 
node, and the following may be written: 


where A is the average boundary surface area associated 
with the mesh point. In Eq. (39), the reason for division 
by Dj is to reduce the coefficients of strains in the 
stress boundary equations to the same order of magni- 
tude as those of Eq. (33). The procedure described here 
for a three-dimensional solid may be readily extended to 
other types of structures with the help of previous 
paragraphs. 
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III. COMMON Variables and COMMON Blocks 

of the Program 


The memory organization in each of the four links of 
ELAS is illustrated in Fig. III-l. Table III-l lists the 
blocks of COMMON sequentially, and gives a short 
description of each block. In the table, the variable- 
address blocks are listed with increasing COMMON ad- 
dresses. It should be noted that the variable-address 
blocks in COMMON are packed in a string, one after the 
other, without any waste of core locations. Such blocks 
may be properly located by means of pointers , which are 
also in COMMON. A pointer is a word whose content 
is one less than the COMMON address of the first word 


of the associated COMMON block. The constituents of 
Block Group 1 are listed in Table VI-3 (Vol. I), in the 
order in which they appear in COMMON. These con- 
stituents are alphabetically ordered with their symbolic 
names in Table III-2. In Table III-3, the meanings of 
entries of important vectors, especially those defined by 
the pointers, are given. The additional COMMON vari- 
ables of Link 4 are listed alphabetically in Table III-4, 
and with increasing COMMON addresses in Table III-5. 
Table III-5 also contains a short description of these 
variables. 
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Table IIM. Sequence and descriptions of COMMON blocks' 


Legend 

** variable address blocks with unique pointers are 
listed with increasing COMMON addresses 

t pointer is a word whose content is the COMMON 
address of the word immediately before the block 
(see Table ill-2 for COMMON addresses of pointers) 

i see Table HI-2 for the meaning of symbols and their 
COMMON addresses 

§ G generate for next link 
S generate for immediate use 
1 inherit from previous link and preserve 
M inherit from previous link and modify 

A see Table VI-3 (Vol. 1) and Table 111-2 for details 
A see Tables HI-4 and II 1-5 for details for Link 4 
□ see Table IV-3 (Vol. 1) 

■ see Table VI-2 (Vol. 1) 

O Li = 2, 9, or 21, depending upon whether ITYPE is 0, 
1 , or 2, respectively 

•■lea — 1, 1, or 3, depending upon whether ITYPE is 0, 
1 , or 2, respectively 

§UOI4DJ9dO 

P >!“!! 

i/> 

— co — c/j co to * — 


- 




§uoi|Djado 

€ 

- - O “ 


% 

- 

3 

to 

§uoijDjedo 

Z 4“!! 

-W"W«WtO- 


Z 

- 

O 


§uo||Djedo 

i 4 u n 

O O O 

0000000000000 000000000000000 

O 

O 



Jt -C 

S m 
to g 

17 

10 

52 

28 

24 

69 

144 

6 

CO <3 

9 OOOOO < 1 ? 

co oooooQQ^ _ ° 9 

O i— f- i— b- >- i— i— i— i- i— Z Z Op •** ol7WWw<<??5ZZZZ 

IND1 

ISUM + 1 

IORD1 

ISUM 

Common address 
of first word 
of block 

1 

17 

27 

79 

107 

131 

200 

344 

351 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Variable 

Pointer 
of the 
blockt 


IIS 

J1 

J2 

J3 

J4 

J5 

J6 

J 7 

JS 

J9 

J10 

IBB 

IBO 

IID 

IlA 

(PR 

ITE 

IDT 

IDY 

IDZ 

ICAR 

ICIX 

I'CIY 

IC1Z 

1C FI 

IXX 

IYY 

IZZ 

IIC 

IL 

Ul 

Q 

9 


MAX 

Block name 

Important constants A 
General descriptors |A 
Important constants and pointers |A 
General descriptors llA 
Subelement load vectorA 
General descriptors III A 
General descriptors IVA 
important constants and pointers ilA 

Subelement stiffness matrix 
J1W vector □ 

J2W vector □ 

J3W vector □ 

J4W vectorD 
J5W vectorD 
J6W vectorD 
J7W vectorD 
J8W vectorD 
J9W vectorD 
J10W vectorD 
IBB vector ■ 

IBO vector ■ 

Material constants vector 
Thermal expansion constants vector 
Pressures vector 
Thicknesses vector 

Uniform temperature increases vector 

Temperature gradients (/-direction) vector 

Temperature gradients (z-direction) vector 

Cross-sectional areas vector 

Torsion constants vector 

Moments of inertia (about y-axi$) vector 

Moments of inertia (about z-axi$) vector 

Angle types vector 

X-coordi nates of nodes vector 

Y-coordi nates of nodes vector 

Z-coordi nates of nodes vector 

C-constants vector® 

Deflections (or loads) vector 

Diagonal-element counts vector 

Reduced stiffness matrix vector 

MAX vector 

sdnojfi 

313013 

- 

<s 

CO 

-r 

_ 

«o 

*o 
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Legend 

** variable address blocks with unique pointers are 
listed with increasing COMMON addresses 

1 pointer is a word whose content is the COMMON 
address of the word immediately before the block 
(see Table ill-2 for COMMON addresses of pointers) 

1 see Table Ilf-2 for the meaning of symbols and their 
COMMON addresses 

§ G generate for next link 
S generate for immediate use 
1 inherit from previous link and preserve 
M inherit from previous link and modify 
♦ see Tables IN-4 and 111-5 
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- 



to 
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e >i“n 
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§uoi|pjedo 
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tO (O to to (O to to 



Block 

length* 

Q 

z 

5000 

5000 

8100 

540 

540 

540 

540 

5994 


Common address 
of first word 
of biock 

Variable 

09000 

09000 

09000 

17100 

17640 

17640 

18180 

00060 


Pointer 
of the 
blockt 

s— 

to 




Block name 

Residual forces vector 

DUMMY vector 
IDUM vector 
ABIN matrix 
ISIR vector 
IMAX vector 
ISIRT vector 
IMIN vector 

BB matrix 

♦ 

V. 

o 

u 

O 

> 

ti- 

ll. 

sdnojB 
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Table III-2. Alphabetical listing of the constituents of COMMON block group l a 


Symbol 

Location in 
COMMON 

Brief description 


Symbol 

Location in 
COMMON 

Brief description 

AA 

l— * * * 

Name of whole COMMON block for floating- 


IELT 

28 

Element type number 



point references 


IERR 

79 

Error indicator 

ACEL 

39 

Body force per unit volume 


IGEM 

78 

Indicator for structures inscribed in 

AL1 

83 

Thermal expansion coefficient of an element 




(Z “ 0)-plane 



in first material axis direction 


IH 

10 ! 

Maximum number of vertices 

AL2 

84 

Thermal expansion coefficient of an element 
In second material axis direction 


1 1 A 

62 

Pointer for thermal expansion coefficient array 

AL3 

85 

Thermal expansion coefficient of an element 


IIC 

74 

Pointer for dbc unit constants array 



in third material axis direction 


(ID 

61 

Pointer for material constants array 

CONS 

45 

Constant for element load vector 


IIS 

77 

Pointer for subelement stiffness matrix 

DG 

82 

Temperature gradient for an element in 
direction y (or z) 


IMAT 

7 

(IIS = 350) 

Number of material types 

DGY 

332 

Temperature gradient along local y-axis for 


IMES 

326 

Indicator for mesh topology input 



an element 


IMET 

31 

Material type number 

DGZ 

331 

Temperature gradient along local z-axis for 


IMF! 

15 

Number of angle types 



an element 


IMMX 

12 

Number of torsion constants types 

DT 

81 

Value of temperature change for an element 


1MMY 

13 

Number of y-moment of inertia types 

D21 

86“ 106 

Material constants for an element 


IMMZ 

14 

Number of z-moment of inertia types 

G1 

47 

First direction cosine of acceleration vector 


IMS 

34 

Number of vertices of current element 

G2 

48 

Second direction cosine of acceleration vector 


IN 

1 

Total number of nodal points 

G3 

49 

Third direction cosine of acceleration vector 


1ND 

33 

IND - IDEG * IN 

IA 


Name of whole COMMON block for fixed- 
point references 


INP 

42 

Indicator for output level 

IARE 

16 

Number of cross-sectional area types ( 


1NX 

9 

Number of last link to be executed 

IBB 

59 

Pointer for IBB array 


IORD 

37 

Number of words allocated for the reduced 
stiffness matrix 

IBN 

2 

Total number of dbc input units 


IOD1 

38 

IORD1 = IORD + 1 

IBO 

60 

Pointer for IBO array 

| 


IP 

4 

Total number of nonzero concentrated load 

IBUN 

327 

Indicator for boundary conditions input 




components 

ICAR 

66 

Pointer for cross-sectional areas array 

| 

IPBG 

43 

Integer constant for element load vector 

ICFI 

70 

Pointer for angles array 


IPEN 

44 

Integer constant for element load vector 

ICIX 

67 

Pointer for torsional constants array 


IPIR 

329 

Indicator for local coordinate axes selection 

1CIY 

68 

Pointer for y-moments of inertia array 


1 PR 

333 

Pointer for pressure array 

1CIZ 

69 

Pointer for z-moments of inertia array 


IPRS 

5 

Number of pressure types 

ICOR 

328 

Indicator for coordinates input 


ISDT 

348 

Number of temperature change types 

3DEF 

75 

Pointer for unknown deflections array 


ISDY 

347 

Number of temperature gradients along local 



(initially loads array) 




y-axis 

1DEG 

8 

Degrees of freedom at a node 


ISDZ 

346 

Number of temperature gradients along local 
z-axis 

IDS 

36 

Order of the subelement stiffness matrix 


ISHUF 

35 

Relabeling indicator 

IDT 

63 

Pointer for temperature changes array 


1ST 

76 

Pointer for reduced stiffness matrix of the 

IDY 

64 

Pointer for temperature gradients array 




whole structure 



(y-direction) 


ISTR 

27 

Indicator for plane-strain case 

IDZ 

334 

Pointer for temperature gradients array 


JSUM 

32 

Number of equations in the reduced system 



(z-direction) 


IT 

3 

Total number of elements 


a See Table HI-1 for sequence and descriptions of COMMON blocks. Table ill-2 is a reordering of Table VI-3 (Vol. I), in both of which the “General Descriptors IV" sec- 
tion of the block is partly excluded . 
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Table 111-2 (contd) 


Symbol 

Location in 
COMMON 

Brief description 


Symbol 

Location in 
COMMON 

Brief description 

ITAP 

41 

Chain program tape number 


J4 

53 

Pointer for J4W array 

ITAS 

335 

Scratch tape number 


J5 

54 

Pointer for J5W array 

ITE 

65 

Pointer for thicknesses array 


J6 

55 

Pointer for J6W array 

ITEM 

29 

Temperature change type number 


J 7 

56 

Pointer for J7W array 

me 

30 

Thickness type number 


J8 

57 

Pointer for J8W array 

ITYPE 

6 

Indicator for materia! type 


J9 

345 

Pointer for J9W array 

IU 

46 

Pointer for diagonal-element count vector 


J10 

344 

Pointer for J10W array 

IXX 

71 

Pointer for X-coordinates array 


M 

25 

Label of current element 

IYY 

72 

Pointer for Y-coordinates array 


Ni 

17-24 

Labels of vertices of an element 

IZZ 

73 

Pointer for Z-coordinates array 


NT1C 

349 

Number of thickness types 

18 

11 

Maximum number of words to describe an 


P 

107-130 

Load vector of a subelement 



element 


PRES 

330 

Pressure value for an element 

JARE 

340 

Type number of cross-sectional area 


S 

351--** 

Subelement stiffness matrix 

JMFI 

15 

Number of angle types 


TE 

80 

Value of thickness for an element 

JMMX 

339 

Type number of the torsional constant about 
local x-axis 


UV 

131-154 

Deflections due to temperature changes for 
an element 

JMMY 

338 

Type number of the sectional moment of inertia 


X 

155-162 

Overall X-coordinates for vertices of an element 



about local y-axis 


XD 

179-185 

X-coordinates of vertices, other than the first 

JMMZ 

337 

Type number of the sectional moment of inertia 




vertex, of an element relative to the first 



about local z-ax is 




vertex 

JPRS 

343 

Type number of pressure 


Y 

163-170 

Overall Y-coordinates of vertices of an element 

JSDY 

342 

Type number of temperature gradient along 
local y-axis 


YD 

186-192 

Y-coordinates of vertices, other than the first 
vertex, of an element relative to the first 

JSDZ 

341 

Type number of temperature gradient along 
local z-axis 


Z 

171-178 

vertex 

Overall Z-coordinates of vertices of an element 

J1 

50 

Pointer for J1 W array 


ZD 

193-199 

Z-coordinates of vertices, other than the first 
vertex, of an element relative to the first 

J2 

51 

Pointer for J2W array 




vertex 

J3 

52 

Pointer for J3W array 


ZGEM 

40 

Floating-point equivalent of 1GEM 
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Table 111-3, Meanings of the entries of important vectors 


Vector in COMMON 

Meaning of rth entry of the vector (ail divisions are in integer arithmetic sense) 

AA vector 

The rth component of the total COMMON vector in floating-point mode 

D21 vector 

The rth component of a row-listed upper 6 X 6 material matrix (see Fig. Ill-2b, Vof. 1), if if exists 

IA vector 

The rth component of the total COMMON vector in fixed-point mode 

IBB-pointer-related vector 

IBB value of Jth degree of freedom direction at ith node (user's label); i = 1 + (r — 1)/IDEG, 
J = r — (i — 1)* IDEG (see Table VI-2, Vol. 1) 

IBO-pointer-related vector 

IBO value of Jth degree of freedom direction at ith node (user's label); i = 1 + (r — 1J/IDEG, 
J “ r — {/ — 1)* IDEG (see Table VI-2, Vol. 1) 

ICAR-pointer-related vector 

Value of rth-type cross-sectional area, if it exists 

ICFI-pointer-related vector 

Value of rth-type angle defining principal axes of cross section, if it exists 

ICIX-pointer-related vector 

Value of rth-type torsional constant, if it exists 

1 Cl Y-poi nter-related vector 

Value of rth-type y-moment of inertia, if it exists 

1C IZ-po inter-related vector 

Value of rth-type z-moroent of inertia, if it exists 

IDEF-poi nter-related vectors 

(1) Value of prescribed concentrated load in Jth degree of freedom direction at node i (user's 
label); i - 1 -h (r — 1J/IDEG, J = / — (/ — 1)* IDEG 

(2) Value of rth component of reduced load vector (the right-hand-side vector in Fig. 11-1, Vol. 1) 

(3) Value of rth component of reduced deflection vector ({d} vector in Fig. JI-1, Vol. 1) 

(4) Value of deflection at the Jth degree of freedom direction at node / (user's label); i “ 1 + 
(r ~ 1)/IDEG, J - r (i - 1)* IDEG 

For (2) and (3) the node i (user's label) and direction J associated with the rth entry may be obtained 
as follows; let r" be the entry number of the word, in IBB-pointer-related vector, where the abso- 
lute value is r and r"th entry of IBO-pointer-related vector is — 1. Then i — 1 + (r" — 1)/IDEG, 
and J •== r" - (i - 1)* IDEG 

IDT -pointer-related vector 

Value of rth-type temperature increase, if it exists 

IDY-poi nter-related vector 

Value of rth-type temperature gradient in /-direction, if if exists 

1 DZ-poi nter-related vector 

Value of rth-type temperature gradient in z-d irection, if it exists 

HA-pointer-related vector 

Value of thermal expansion coefficient in the Jth material axes direction in element i; i - I + 
(r ~ 1)/k 2 , J — r — ks* (i — 1), where k 2 is 1, 2, or 3, depending upon whether ITYPE is 0, 1, 
or 2, respectively 

IIC -pointer-related vector 

Value C of Jth degree of freedom direction at ith node (user's label); / = 1 + (r — 1J/IDEG, 
J = r - (i ~ 1)* IDEG (see Table Vl-2, Vol. 1) 

11 D-poi nter-related vector 

Value of Jth material constant of material type i; 1=1 + (r — 1)/ki, J = r — ki* (i — 1), where 
ki is 2, 9, or 21, depending upon whether ITYPE is 0, 1, or 2, respectively 

HS-pointer-reJoted vector 

Element kmn of the free-free subelement stiffness matrix, m = 1 4- (r — 1J/IDS', n == r — IDS'* 
(m — 1); m corresponds to m'th degree of freedom direction (m' = 1 + (m — l)/IMS f ) at 
vertex m" (m" = m — IMS'* (m — 1); n corresponds to n'th degree of freedom direction 
(n # ■= 1 + (n — 1)/1MS') at vertex n" (n^ = n — IDS 1 '* (n' — I); IMS 7 is the number of vertices 
of subelement, and IDS' = IMS'* IDEG 

1 PR-pointer-related vector 

Value of rth-type pressure, if it exists 

IST-pointer-related vector 

(1) Element Km n of the stiffness matrix of the supported structure. To find mth direction, enter IBB- 
pointer-related vector with the entry number / of the word, in lU-pointer-refated vector, which 
is closest to, but not greater than r. let r" be the entry number of the word, in IBB-pointer- 
related vector, whose absolute value is r and the r"th entry in IBO-pointer-reJafed vector is 
— 1; mth direction corresponds to m"th degree of freedom direction at node m' (user's label); 
m' = 1 -j- (r" — 1J/IDEG, m" = r" — (m' — 1)* IDEG. To find nth direction, determine s' by 
adding to r the difference between r and the r'th entry of lU-pointer-related vector. Let s" be 
the entry number of the word in IBB-pointer-related vector, whose absolute value is s' and the 
s"th entry in IBO-pointer-related vector is " — 1 ; nth direction corresponds to n"th degree of free- 
dom direction at node n (user's label); n' = 1 4* (s" — 1)/1DEG, n" = s" — [n — 1)* IDEG 

(2) Value of residual force acting at node i in direction J, where i ~ 1 4* (r — 1J/IDEG, j = r — 
(i - 1)* IDEG 
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Table 111-3 (contd) 


Vector in COMMON 

Meaning of rth entry of the vector (all divisions are in integer arithmetic sense) 

ITE-pointer-related vector 

Value of rth-type thickness, if it exists 

lU-pointer-related vector 

Entry number in IST-pointer-related vector of rth diagonal element of the reduced stiffness matrix 

1 XX- pointer-related vector 

X-coordinate of node r (user's label) 

lYY-pointer-related vector 

Y-coordinate of node r (user's label) 

IZZ-poi nter-related vector 

Z-coordinate of node r (user's label), if it exists 

J 1 -poi nter-related vector 

J1W value of rth element (see Table IV-3, Vol. 1) 

J2-pointer-related vector 

J2W value of rth element (see Table IV-3, Vol. 1) 

J3-poi nter-related vector 

J3W value of rth element (see Table IV-3, Vol. 1) 

J4-poi nter-related vector 

J4W value of rth element (see Table IV-3, Vol. 1) 

J5-pol nter-related vector 

J5W value of rth element (see Table IV-3, Vol. 1) 

J6-poi nter-related vector 

J6W value of rth element (see Table IV-3, Vol. 1), if it exists 

J7-pointer-reIated vector 

J7W value of rth element (see Table IV-3, Vol. 1), if it exists 

J8-pointer-related vector 

J8W value of rth element (see Table IV-3, Vol. 1), if it exists 

J9-pointer-related vector 

J9W value of rth element (see Table IV-3, Vol. 1), if it exists 

J 10-pointer-related vector 

J10W value of rth element (see Table IV-3, Vol. 1), if it exists 

N vector 

The label (user's) of the rth vertex of an element 

MAX-poi nter-related vector 

Number of nonzero entries above rth diagonal element of the decomposed reduced stiffness matrix 
(see Fig. 1 1-1, Vol. 1) 

P vector 

Element load acting in direction J of ith vertex of a subelement; J = 1 + (r — 1)/1MS 7 , i — r — 
(J — 1)* IMS' (IMS' — number of vertices of the subelemenf) 

S matrix 

See HS-pointer-related vector 

UV vector 

Deflection in direction J of ith vertex of a subelement subjected to temperature change in local coor- 
dinates; J — 1 + (r — 1)/IM$', i = r — (J — 1)* IMS 7 (IMS 7 — number of vertices of subelement) 

X vector 

X-coordinate of rth vertex of an element 

XD vector 

X-coordinate, relative to the first vertex, of (r + l)st vertex of an element 

Y vector 

Y-coordinate of rth vertex of an element 

YD vector 

Y-coordinate, relative to first vertex, of (r + l)st vertex of an element 

Z vector 

Z- coordinate of the rth vertex of an element 

ZD vector 

Z-coordinate, relative to the first vertex, of (r + 1) st vertex of an element 
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Table 111-4. Alphabetical list of additional COMMON variables for Link 4 a 


Symbol 

COMMON location 


Symbol 

COMMON location 


Symbol 

COMMON location 


Symbol 

COMMON location 

A 

14785— 15415 


ICAS 

212 


IWG 

14660-14749 


NU 

292-294 

ANGLE 

211 


ICLA 

206 


JM1 

293 


QF 

253-258 

ARE 

205 


ICLAS 

274-277 


JP1 

292 


QN 

247-252 

AST 

203 


ICN 

201 


JS1 

294 


RED 

265-270 

B 

15416-15479 


ICOL 

295 


LM 

202 


RES 

259-264 

BAS 

271-273 


ICON 

210 


MAC 

U400-14659 


SIR 

223-225 

SIR 

220-222 


IDR 

297 


MSET 

15596-15695 


SR 

235-240 

BST 

217 


IE 

213 













MB 

215 


W 

15696-15704 

C 

15480-15495 


IM 

208 







DD 

14750-14785 


IMEL 

207 


NB 

214 


XF 

244-246 

DIN 

226-234 


INBON 

204 


NBAN 

278-287 | 

i 

XN 

241-243 

ETA 

229-231 

j 

IONE 

200 


NEl 

14000-14399 

j 

XII 

226-228 

FF 

14000-15704 


IR1G 

296 


NES 

295-297 


ZTA 

232-234 

1C 

209 


IROT 

216 


NSET 

15496-15595 




a See Table II 1-5 for meanings of variables. 


Table 111-5. List of additional COMMON variables for Link 4 a 


Location in COMMON 

Symbol 

Brief description 

1-199 


This portion of COMMON is as in Table VI-3 of Vol. 1 

200 

IONE 

Total number of one-dimensional elements in the structure 

201 

ICN 

Label of current mesh point (ICN varies from 1 to IN) 

202 

LM 

Total number of non-one-dimensional elements meeting at mesh point ICN 

203 

AST 

Indicator containing * or BCD blank, depending upon whether mesh point ICN is on boundary or not, 
respectively 

204 

INBON 

Indicator containing 1 or 0, depending upon whether mesh point ICN is on boundary or not, respectively 

205 

ARE 

Average boundary surface area for mesh point ICN, if it is on boundary 

206 


Total number of class** types for elements of material type group IM at mesh point ICN 

207 

1 ■ 

Total number of material types at mesh point ICN 

208 


Current material type group number (IM varies from 1 fo IMEL) 

209 

1 

Current class* type group number (1C varies from 1 to ICLA) 

210 

■ 

Sequence number of current strain- deflection equation at mesh point ICN for material group IM and for 
das s b group 1C 

211 

ANGLE 

Angle between XII local axis and the 1—2 line of the lowest labeled shell element attached to mesh 
point ICN 

212 

ICAS 

Cla$s b type number of ICth cla$s b group of IMth material group at mesh point ICN 

213 

IE 

Number of mesh elements (of class b group 1C of material group J M) plus 1 at mesh point ICN 

214 

NB 

Total number of mesh points in node set at mesh point ICN 

215 

MB 

Number of boundary points attached to mesh point ICN 

216 

IROT 

Indicator containing 0 or 1, depending upon whether local axes at mesh point ICN are parallel to overall 
axes or not, respectively 

217 

BST 

Indicator containing BCD blank or **, depending upon whether local axes at mesh point ICN are parallel 
fo overall or not, respectively 

220-222 

BIR 

Direction cosines of outer unit normal vector at mesh point ICN, if it is on boundary 

a Th|s fable is not applicable fo subroutine DIM! of Link 4. 
b CIass types are those of Table Vl-6, Vol. 1. 
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Table !ll-5 (contd) 


Location in COMMON 

Symbol 

Brief description 

223-225 

SIR 

Vector heading towards structure at mesh point ICN, if it is on boundary 

226-234 

DIN 

Direction cosines of local axes In overall coordinate system at mesh point ICN (the columns of DIN are 
named as XII, ETA, and ZTA) 

235-240 

SR 

Independent components of stress tensor for ICth class* group of IMth material group at mesh point ICN 

241-243 

XN 

Overall coordinates of mesh point ICN 

244-246 

XF 

Overall coordinates of llth vertex of iLth element of ICth class* group of IMth material group at mesh 
point ICN 

247-252 

QN 

Deflection components in overall coordinates of mesh point ICN 

253-258 

QF 

Deflection components in overall coordinates of mesh point whose overall coordinates are in XF 

259-264 

RES 

Residual forces* in overall coordinates at mesh point ICN, if on boundary 

265-270 

RED 

Relative deflections (in overall coordinates) of mesh point related with XF vector with respect to mesh 
point ICN 

271-273 

BAS 

Direction cosines of 1-2 line of the lowest labeled element of class* group 1C of material group IM at 
mesh point ICN 

274-277 

ICLAS 

Number of class* groups in each material group (maximum 4) of mesh point ICN 

278-287 

NBAN 

List of labels of boundary mesh points attached to mesh point ICN 

292-294 

NU 

Vector containing the sequence numbers of the vertices after (JPl), before (JM1), and above (JS1) mesh 
point ICN in the 1th mesh element of the node set (with Table 111-5, vol. 1) 

292 

JPl 

See NU (1) 

293 

JM1 

See NU (2) 

294 

JS1 

See NU (3) 

295-297 

NES 

Vector containing number of independent strain components (ICOL), number of right-hand sides (IRIG) 
and indicator of right-hand-side arrangement (IDR) (IDR = 0 means lineal strains first, IDR — 1 means 
rotational strains first) for current ICN/IM/IC 

295 

ICOL 

See NES (1) 

296 

IRIG 

See NES (2) 

297 

IDR 

See NES (3) 

329—349 


This portion of COMMON is as shown in Table VI-3, Vol. 1. See also Fig. 111-3, Link 4 

349-13999 


See Fig. ill-1, Link 4 

14000-15704 

FF 

Vector containing information for stress conputation at mesh point ICN 

14000-14399 

NEL 

Element set information of mesh point ICN (see Table VI-7, Vol. 1) 

14400-14659 

MAC 

Table for classes and material of element set at mesh point ICN (see Table VI-7, Vol. 1) 

14660-14749 

IWG 

Vector of weights of strain-deflection equations for current ICN/IM/IC 

14750-14785 

DD 

i 

Material matrix for current ICN/IM/IC 

14786-15415 

A 

Augmented matrix of strain-deflection equations for current ICN/IM/IC 

15416-15479 

B 

Coefficient matrix (or its inverse) of the least-squares equations for strain for current ICN/IM/IC 

15480-15495 

C 

Right-hand-side vectors) of the least-squares equations for strains for current ICN/IM/IC 

15496-15595 

NSET 

List of labels of mesh points on the boundary and attached to mesh point ICN 

15596-15695 

MSET 

Auxiliary array for NSET 

15696-15704 

W 

Direction cosines of new material axes in the old for current ICN/IM/IC 


c Residuai forces are those listed by Output Item 20 (see Sect. Vl-D and Vl-E, Vol. I). 
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